Dynamics of an Inelastic Gravitational Billiard with Rotation 
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The seminal physical model for investigating formulations of nonlinear dynamics is the billiard. 
Gravitational billiards provide an experimentally accessible arena for their investigation. We present 
a mathematical model that captures the essential dynamics required for describing the motion of a 
realistic billiard for arbitrary boundaries, where we include rotational effects and additional forms 
of energy dissipation. Simulations of the model are applied to parabolic, wedge and hyperbolic 
billiards that are driven sinusoidally. The simulations demonstrate that the parabola has stable, 
periodic motion, while the wedge and hyperbola (at high driving frequencies) appear chaotic. The 
hyperbola, at low driving frequencies, behaves similarly to the parabola; i.e., has regular motion. 
Direct comparisons are made between the model's predictions and previously published experimental 
data. The representation of the coefficient of restitution employed in the model resulted in good 
agreement with the experimental data for all boundary shapes investigated. It is shown that the 
data can be successfully modeled with a simple set of parameters without an assumption of exotic 
energy dependence. 

PACS numbers: 05.10.-a;05.45.a;05.45.Ac;05.45.Pq 



I. INTRODUCTION 



The seminal physical model for investigating formula- 
tions of nonlinear dynamics is the billiard. Tradition- 
ally this is a field free, classical, system where a parti- 
cle experiences elastic collisions with a rigid boundary. 
Depending on the boundary shape, the ensuing motion 
can be stable or chaotic [11, Q . Quantum mechanical bil- 
liards have also been investigated and many of their es- 
sential features can be effectively represented with clas- 
sical waves on membranes or cavities . A limitation of 
the classical billiard is the difficulty in reproducing the 
system in the laboratory as a consequence of dissipation 
and the earth's ubiquitous gravitational field. In con- 
trast, gravitational billiards provide an experimentally 
accessible arena for testing formulations of nonlinear dy- 
namics. 

One and two-dimensional Hamiltonian versions of 
gravitational billiards have long provided easily visual- 
ized systems that exhibit a wide range of stable and 
chaotic behavior j^-Q- The system consists of a par- 
ticle undergoing elastic collisions within a rigid bound- 
ary, where the particle follows a ballistic trajectory un- 
der the influence of a constant gravitational field between 
collisions. When the boundary is periodically driven, 
Fermi acceleration may result [9] , establishing a connec- 
tion with cosmic ray physics and cosmology. Billiards 
that bounce vertically on a level, oscillating surface have 
been used to model the impact process for numerous en- 
gineering applications including moving parts in machin- 
ery, impact dampers, fluid induced vibration in tubes 



and moored ships driven by steady waves 0, [T^ [TT| . Re- 
cent interest in dynamics has been focused on dissipa- 
tive systems such as granular media. While inelasticity 
in these systems is usually represented by a collisional 
restitution coefficient, it has been observed that the in- 
clusion of rotational friction induces qualitative changes 
in behavior [12] that cannot be explained by other means. 
Similarly, in billiard experiments [13| , when friction is left 
out of the theoretical formulation, it appears that one is 
forced to make unphysical assumptions about the source 
of energy loss to approximately replicate the experimen- 
tal data 
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This paper considers the more realistic situation of an 
inelastic, rotating, gravitational billiard in which there 
are retarding forces due to air resistance and friction. In 
this case the motion is not conservative, and the billiard 
is no longer a particle, but a sphere of finite size. We 
present a mathematical model that captures the relevant 
dynamics required for describing the motion of this "real 
world" billiard for arbitrary boundaries. The model is 
applied to parabolic, wedge and hyperbolic billiards that 
are driven sinusoidally. Direct comparisons are made be- 
tween the model's predictions and experimental data pre- 
viously collected [lj|- Although several studies have in- 
vestigated the effect of variable elasticity in relation to 
the gravitational billiard, this study is the first to incor- 
porate rotational effects and additional forms of energy 
dissipation. 

The ergodic properties of Hamiltonian gravitational 
billiards are well studied [l-d, [l^. It has been shown 
that the parabolic billiard is completely integrable hav- 
ing stable, periodic orbits Q. Studies of the wedge bil- 
liard demonstrated that the billiard's behavior depends 
on the vertex angle, defined as 20 0. For < < 45°, 
the phase space contains coexisting stable and chaotic be- 
havior. For Q — 45° the motion is completely integrable, 
while for B > 45° the motion is chaotic. Wojtkowski 
refers to this geometry as "fat billiards" and has rigor- 
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oiisly proven that they have a single, ergodic component 
[Tsj . These resuhs have also been confirmed through ex- 
periments for an optical billiard with ultra cold atoms 
|16| . It has been demonstrated with numerical simulation 
that the motion of a hyperbolic billiard exhibits charac- 
teristics of the parabolic billiard for low energy, where 
the motion is near the origin, and for the wedge billiard 
at high energy, where the motion is mostly concentrated 
at its asymptotic limits [Ij. 

Feldt and OlafsenfT!] have experimentally studied a 
real inelastic billiard for a variety of boundaries. One 
experiment consisted of a steel ball moving within a 
closed reflective aluminum boundary shaped either as 
a parabola, wedge or hyperbola. The container was 
driven in the horizontal direction to compensate for en- 
ergy losses resulting from collisions. Imaging software 
determined the ball's position and velocity at the colli- 
sion points. The study results indicated regular motion 
for the parabola and chaotic motion for the wedge; the 
motion for the hyperbola was found to be frequency de- 
pendent, sharing characteristics of the parabola at low- 
driving frequencies and the wedge for higher-driving fre- 
quencies. 

In this work direct comparisons are made between sim- 
ulations of the model system and the experimental data 
of Feldt and Olafsen. To date, Gdrski and Srokowski[l^ 
are the only investigators known to have theoretically 
studied the experiments conducted by Feldt and Olaf- 
sen. In their model they consider an inelastic, gravita- 
tional billiard for the case of no friction (or rotation) and 
no drag. In order to replicate the main features of the 
experiments, it was necessary to resort to a surprising, 
unconventional representation of the restitution coeffi- 
cient energy dependence. 

This paper begins with a discussion on Kane's equa- 
tions and the impact theory used for describing a col- 
lision between a billiard and a moving boundary. The 
Appendix develops the equations governing this collision 
process in detail. Sections III and IV describe the trajec- 
tory model for tracking the billiard's motion after each 
bounce and the procedure for detecting collisions. Sec- 
tion V explains how the coefficients of restitution and 
friction are determined for numerical simulations. Sec- 
tion VI presents simulations comparing the numerical re- 
sults to previous experiments. Conclusions are presented 
in Section VII. 



II. KANE'S EQUATIONS AND THE IMPACT 
THEORY 

There are several competing models for describing a 
collision of a billiard with a boundary, and each model 
has its own merit depending on the intended applica- 
tion. First, there are analytical models that determine 
the billiard's velocity post collision in terms of the pre- 
impact velocities. The formulation is based on Newton's 
law of motion and Coulomb's law of friction, and requires 



knowledge (a priori) of the coefficients of restitution and 
friction. Accepted models in this category include works 
by Waltonfm and Kane and Levinson[l8|. Second, there 
are impact models based on the field of continuum me- 
chanics, which consider collisions of elastic, viscoelastic 
or plastic objects. Models in this group are based on the 
Hertzian contact theory and its offshoots fioj. 

In this study we employ a modified version of the im- 
pact theory set forth by Kane and Levinson for collisions 
between the billiard and boundary where we account 
for the effects of a moving boundary. The original the- 
ory provides a direct method for computing the billiard's 
generalized speeds post collision considering fixed bound- 
aries. The theory is based on Kane's equations which 
utilize partial velocities and generalized forces for deriv- 
ing equations of motion. The equations are also known 
as La gran ge's form of D'Alembert's principle, and refer- 
ences [T8!| and f20| provide a thorough treatment on the 
subject. 

Kane and Levinson make the following assumptions in 
their model: first, the contact area between the objects 
is a single point through which all forces are exerted. 
Second, the total collision impulse is represented by the 
integral of the forces over the entire collision time. Third, 
the coefficients of restitution, static friction and kinetic 
friction are constants to be determined experimentally. 
The theory initially assumes no slipping at the contact 
point between the sphere and boundary. A set of values 
for the generalized speeds are generated, and are valid if 
and only if the no-slip condition is satisfied. If the no- 
slip condition is violated, then a new set of values for the 
generalized speeds are developed under the assumption 
of slipping. See reference [l^ for a detailed derivation of 
the theory. 

We extend Kane and Levinson 's impact theory to in- 
clude collisions on moving boundaries. For completeness, 
we list the equations of motion for this system in the 
Appendix. Next, we introduce a trajectory model that 
tracks the billiard's motion between bounces, taking into 
account dissipative and aerodynamic forces. By employ- 
ing the impact and trajectory models, we then construct 
an efficient set of algorithms describing the billiard's mo- 
tion in order to perform numerical simulations. 



III. TRAJECTORY MODEL AND THE 
REINITIALIZATION OF THE GENERALIZED 
SPEEDS 

Between collisions we make use of a trajectory model 
that numerically tracks the billiard's motion after each 
bounce. The model is used to reinitialize the generalized 
speeds at the point of initial contact with the boundary. 
As a starting point, we define the equations of motion 
that governs the billiard's movement while airborne. 

For a billiard (or sphere) moving through air, its mo- 
tion is affected by gravity, air resistance (drag) and ad- 
ditional aerodynamic forces due to its spinning motion 
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(the Magnus effect). In this study the Magnus Force is 
neglected since its overall effect is small. The force of 
gravity acting on the billiard is defined as 



(1) 



where m is the billiard's mass and g is the acceleration 
due to gravity. Refer to the Appendix for a definition of 
the coordinate system. 

The drag force exerted on the billiard is a function 
of the billiard's velocity, and acts in the direction op- 
posite to its path. At low velocities the drag force is 
linearly proportional to the billiard's speed, but shows a 
quadratic dependence on speed at higher velocities. Gen- 
erally, the drag force acting on a body is determined by 
experimental measurements, and is often approximated 
by the equation 



Yd = -v(ci 



C2 |v| 



(2) 



where ci and C2 are constants that are dependent on the 
size and shape of the object [2l|. Traditionally, ci and 
C2 are expressed as 



ci = dTTr/b 
C2 = ^pACo 



(3) 



where rj is the dynamic viscosity of air, b is the object's 
radius, p is the atmospheric density, A is the object's 
cross-sectional area and Cd is the drag coefficient. In 
Equation ([3]), ci is the coefficient of the familiar Stokes 
drag force, which is valid at low Reynolds number. The 
drag coefficient, Cd, is a function of the Reynolds num- 
ber. Experimentally, ci and C2 have been measured di- 
rectly, giving the correct dependence on the object's di- 
ameter. For spheres in air, approximate values for ci and 
C2 in SI units are 



ci = 1.55 X 10-^ 
C2 = 0.22D'^ 



(4) 



where D is the sphere's diameter in meters [21|. In this 
study constants ci and C2 are specified by Equation ([4]), 
but note that the constants given by Equation ^ yield 
similar results. The ratio of the quadratic term to the 
linear term of the drag force, i.e.. 



C2V V ^ , , „ 

' ' = 1.4 X 10^ V D 

ClV 



(5) 



determines which type of drag is more significant. If the 
ratio is greater than one, the quadratic term is dominant. 
If the ratio is below 1, the linear term is dominant. If the 
ratio is around one, however, both terms must be taken 
into account. 

Subsequently, the equations of motion for a billiard 
traveling through air is 

F = Fg + (6) 



The billiard and boundary are simulated using a time- 
driven procedure, where the system is advanced in time 
until a collision is detected [IM H^l- Between bound- 
ary encounters the trajectory equations denoted by equa- 
tion ([6]) consist of second-order, nonlinear, coupled dif- 
ferential equations which are solved numerically by using 
a fourth-order Runge-Kutta method. 



IV. COLLISION DETECTION METHOD 

In this section we outline a general procedure for de- 
tecting collisions between the billiard and a boundary 
of arbitrary shape. The procedure locates the minimum 
distance between the objects at each time step, and com- 
pares that distance to a specified tolerance, which for our 
case is the billiard's radius b. If the distance is less than 
or equal to the tolerance, then a collision is reported. 
Otherwise, it is concluded that no collision has occurred. 

A detailed description of the procedure now follows: 
first, write the square of the distance formula be- 
tween the billiard and boundary in terms of the billiard's 
geometric center and the mathematical formula that de- 
scribes the boundary's shape. Second, use the boundary 
formula to eliminate all but one of the variables, thereby 
expressing as a function of a single variable. Third, 
minimize by taking its derivative and setting it equal 
to zero. Fourth, solve for the roots of the resulting equa- 
tion, where valid solutions are restricted to the set of real 
numbers. (Depending on the boundary, the roots may be 
determined by analytical or numerical methods). Fifth, 
locate the minimum distance between the objects by sub- 
stituting the roots into L^. (If , then the billiard 
is impacting the boundary. Otherwise, the objects are 
not colliding). 

If a collision is detected, the collision time and colli- 
sion location are approximated by interpolation methods. 
The procedure for finding the collision time is based on a 
paper by Baraff |25j] . It searches for a configuration where 
the penetration depth between the objects is sufficiently 
close to zero. As a consequence the determination of the 
collision time is transformed into a root-finding problem, 
where the system's state at the collision time is approxi- 
mated by interpolating the derivatives computed by the 
Runge-Kutta method. 

In the following we consider driven parabolic, wedge 
and hyperbolic boundaries defined mathematically (in 
the laboratory frame), respectively, as: 

q2 = .f{qi) = a{qi- Aqif + c (7) 
92 = ./(gi) = 6|(7i-Agi|+c (8) 

f{qi)= sj^l + ^i^^i^^^^-S (9) 

where a — 0.26cm^^, b — 1.85, c = 0.63cm, a = 40.3cm^, 
/3 = 0.08cm~^ and d = AA5cm. These are the values used 
in the experiments of Feldt and Olafsen. The boundaries 
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oscillate horizontally and their position at time t is de- 
fined by 

Aqi {t) = Asinujt (10) 

where A is the amplitude and uj — 27r/ is the oscilla- 
tion angular frequency in rad/sec and / is the oscillation 
frequency in hertz. Experimentally, note that the bound- 
aries are sealed off by an aluminum top and thin pieces 
of Plexiglas on the sides, rendering the system effectively 
two-dimensional. Figure [T] shows the boundary shapes 
and their orientation with respect to gravity and the driv- 
ing direction. Note the boundaries are offset and their 
horizontal tops are omitted for clarity. The billiard's po- 
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FIG. 1: A schematic of the boundaries used in the simulations. 



sition is tracked in time by following its geometric center, 
where its position is defined by 

q = qini + q2n2 + qsua (11) 

For the numerical simulations presented in this study, 
53 = 0. In reality, however, q^ ^ because of slight chat- 
tering in the ria direction. Private communications with 
an author of [1^ reveal that the system noise induced by 
this effect is small. The actual billiard considered in the 
simulations and used in the experiments is a 3.1 mm di- 
ameter steel chrome ball weighing approximately 0.13028 
g- 

Application of the above procedure results in solving 
a second-order algebraic equation for the wedge, a third- 
order algebraic equation for the parabola and a fourth- 
order algebraic equation for the hyperbola. We deter- 
mine the roots of the equations by means of the Newton- 
Raphson method. 



V. THE COEFFICIENTS OF RESTITUTION 
AND FRICTION 

Collisions with the boundary result in energy losses 
stemming from the restitution in the normal direction 
and friction in the transverse direction. Feldt and Olaf- 
sen [l^ suggest a coefficient of restitution of 0.9 between 
the steel billiard and the aluminum boundary, noting 
that its value is velocity dependent. However, as we see 
in the following, the coupling between the normal and 
tangential contact forces during the impact reduces this 
coefficient considerably. Further, it has been observed 
that the coefficient of restitution depends on the incident 
angle at impact |26l] . Friction, due to confining walls, 
also plays a vital role in the experiments of Feldt and 
Olafsen. Studies on granular media show the importance 
of including frictional effects between particles and their 
containers since the particles' velocity distributions are 
affected by these interactions [13, H^. Kane and Levin- 
son's impact theory remains practical when it is supple- 
mented by experimental measurements capturing the co- 
efficients of restitution and friction. Toward this end the 
parabolic billiard is used as a standard test case for es- 
tablishing the coefficient of restitution since experiments 
have shown that it exhibits stable, period-one orbits. The 
experiment described in [l^ resulted in an orbit height 
of approximately 7.5 cm; an apparent value for the coef- 
ficient of restitution is estimated by matching the orbit 
height of the simulation to the experiment. First, note 
that for steel on aluminum, experiments reveal that the 
coefficient of static and kinetic friction is approximately 
0.61 and 0.47 [2^, respectively. If the numerical model 
applies a coefficient of restitution of e = 0.393 along with 
the friction coefficients specified above, then the simula- 
tion approximately replicates the experiment; as a result 
we apply this e value for all boundary shapes considered 
in this paper. Note that if the effects of air resistance are 
omitted, then the coefficient of restitution drops slightly 
to e = 0.392. This result is not surprising considering 
the smallness of the billiard and its relatively short time 
of flight between bounces. Figure [5] shows the trajec- 
tory of the stable orbit and the location of the parabolic 
boundary at the impact points, while Figure [3] reveals the 
evolution of the billiard's trajectory for entering a period- 
one orbit after starting at the origin. If sufficient energy 
is supplied to the system, the billiard's trajectory eventu- 
ally mode-locks into a stable period-one orbit. The orbit 
moves up or down the parabola as the driving frequency 
is increased or decreased, respectively. If insufficient en- 
ergy is given to the system, then the parabola will explore 
multiple trajectories. 

If one ignores friction, the coefficient of restitution re- 
quired to match the orbit height from the experiment 
drops significantly to e = 0.246. If air resistance is 
also neglected, then the coefficient of restitution drops 
to e = 0.245. Table U summarizes the different values of 
the coefficient of restitution considering several dynami- 
cal effects. For the amplitude and frequencies considered 
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in this paper, the coefficient of restitution is dominated 
by frictional effects, while air drag plays a negligible role. 
This demonstrates a greater need to understand the rela- 
tionships between the friction and restitution coefficients. 



TABLE I: The coefficient of restitution required to reach a 
stable, period-1 orbit for a parabolic billiard considering vari- 
ous dynamical effects. X indicates the effect is included; 

indicates the effect is omitted. 



Friction 


Drag 


Coefficient of Restitution 






0.245 




X 


0.246 


X 




0.392 


X 


X 


0.393 



Gorski and Srokowski suggest a coefficient of restitu- 
tion of e = 0.43 for the no friction, no drag case. How- 
ever, their approach for determining the billiard's veloc- 
ity change post collision is unconventional since they ap- 
ply the coefficient of restitution to the complete veloc- 
ity, instead of only to the normal component of velocity. 
Without friction, the parallel component of momentum 
must be conserved. To confirm their assumption, exam- 
ine Equation 1 of their paper, which reproduced here is 

vf =r(vo^-2u(vo^.u)) (12) 

where Vq' and vf are the particle's velocity before and 
after a collision, respectively, u is the velocity normal to 
the boundary and r is the coefficient of restitution 
Now, suppose V is the velocity tangent to the boundary. 
Taking the scalar product of Equation with respect 
to u and v results in the following expressions: 

vf • u = r • u - 2u ■ u (v^ ■ u)] 

= r [v[r • u] [1 - 2] 

= -r[v^u] (13) 

vf • V = r [v^ • V - 2u • V (v^ • v)] 

= ^ [vf? • v] (14) 

By examining Equations (IT^ and ([T^ , it is clear that the 
coefficient of restitution is incorrectly applied to both the 
normal and tangential components of the velocity. As a 
consequence comparison of their results to the experi- 
ments performed by Feldt and Olafsen remain ambigu- 
ous. 

In granular media simulations, a method of preventing 
inelastic collapse of particles is to set the coefficient of 
restitution to its elastic limit of 1 if collisions occur too 
frequently [s^. Additionally, studies have demonstrated 
that the coefficient of restitution approaches a value of 
1 as the normal component of the impact velocity ap- 
proaches 0. As a consequence we apply a coefficient of 



restitution of 1 in our simulations if the relative veloc- 
ity (between the billiard and boundary) in the normal 
direction is sufficiently small that it results in inelastic 
collapse, where a "nearly infinite" number of collisions 
occur in a finite time j30j . In practice, this assumption is 
only applied at the start of simulations (to initiate mo- 
tion) and for brief instances of time during the simulation 
(as explained above). 
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FIG. 2: The trajectory of the stable period-one orbit and the 
location of the parabolic boundary at the impact points. 
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FIG. 3: The evolution of the billiard's trajectory for achieving 
a stable period-one orbit. 



VI. SIMULATIONS 

Each simulation tracks a single trajectory consisting 
of 25,000 billiard-boundary collisions. The billiard is ini- 
tially at rest, but is quickly propelled into the air by 
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the energy transmitted from the boundary to the bil- 
hard. The coUision height 52 and time t between consec- 
utive bounces are extracted from the simulations and are 
shown in Figure HI for multiple boundaries driven at vary- 
ing frequencies. For comparison, the experimental results 
from Feldt and Olafsen [T^] are given in Figure [H The 
successive mappings of the collision heights and times of 
flight show good resemblance to the experimental data 
even though a constant coefficient of restitution is used. 
In reality, the coefficient of restitution is velocity depen- 
dent and its representation is affected by the coupling 
between the normal and tangential contact forces [sl)- 
The numerical model, however, detects additional col- 
lisions (i.e., collisions that occur in rapid succession) not 
reported by the experiments. This is due to the lower 
resolution of the experimental data. As a result the nu- 
merical model observes more collisions at both the small 
and large values of the height (72; thus the time mappings 
have longer time tails associated with shorter times of 
flight for these collisions. 

The plot of the wedge driven at 6.6 Hz reveals that 
the motion appears chaotic as suggested by the experi- 
ments. The billiard is continuously driven to the top of 
the wedge, with most of the collision points laying above 
the 92, n = Q2.n+i Imc in the return map. The time map 
also shows indications of chaotic behavior and a similar 
concentration of points to the experimental data. The 
hyperbola driven at 5.8 Hz resembles the unstable behav- 
ior of the wedge in both position and time, and shows 
a likeness to the experimental results. For lower driv- 
ing frequencies, the hyperbolic billiard is well approxi- 
mated by the parabolic billiard as seen in Figure ID At 
4.5 Hz, the billiard's motion is confined to the regions 
near the hyperbola's vertex, and shows semblance of a 
regular pattern not noted in the experimental data due 
to possible smearing of the data. Patterns are also de- 
tected in the temporal mapping of the hyperbola at this 
driving frequency. Figure |B] displays these patterns in 
both the spatial and temporal mappings for the hyper- 
bolic and parabolic billiards. Note that the parabolic 
billiard driven at 4.5 Hz is not studied experimentally, 
but is used to demonstrate that the hyperbolic billiard 
behaves similarly to the parabolic billiard at low driving 
frequencies. 

Following the experiments the phase space is fur- 
ther investigated by plotting the normalized collision 
height q2/q2,max versus the normalized tangential veloc- 
ity Ui/u4^max post collisiou. The quantities q2,max and 
U4,max are defined as the maximum energy values that 
the billiard can possess at the collision height q2 if all the 
energy were potential or kinetic, respectively. A com- 
pletely stable period-one orbit for a perfectly elastic bil- 
liard is characterized by having zero tangential velocity 
assuming collisions with a symmetric boundary. For the 
parabolic billiard driven at 5.4 Hz, the numerical model 
predicts a small value for the normalized tangential ve- 
locity when the billiard achieves a stable period-one or- 
bit; the mapping is a single point that has the value 



{q2/q2,max,U4/u4,max) = (0.376, ±0.0372). The experi- 
mental data, however, shows that the normalized tangen- 
tial velocity is a thin band about zero, where the range in 
velocity and height is caused by the noise in the system 
and small variations in the coefficient of restitution. Fig- 
ure[7]shows the results for the remaining boundaries. For 
the case of the wedge, the billiard explores much of the 
phase space; the hyperbolic billiard driven at the higher 
frequency exhibits similar behavior to the wedge, but ex- 
amines even more of the phase space. For both shapes 
there are regions that have concentrations of points in the 
phase space not indicated by the experiments. For the 
hyperbolic billiard driven at the lower frequency, regular 
patters develop, which are similar in appearance to the 
parabolic billiard driven at the same frequency. The sim- 
ulations of this paper for the hyperbolic billiard driven 
at 4.5 Hz is comparable to the results reported by Gdrski 
and Srokowski; the similarity exists because for low driv- 
ing frequencies (or for low energy systems) the effects of 
their restitution assumption are mitigated. 

Note that the plots in Figure [7] represents a signifi- 
cant deviation from the results indicated by the experi- 
ments. The difference is qualitative and is not explained 
by the extra collisions reported by the numerical simula- 
tions. Potentially, the source of the difference may lay in 
the finite resolution of the imaging software, resulting in 
measurement uncertainty, where the normalized collision 
height and normalized tangential velocity are calculated 
quantities that depend on the accurate resolution of the 
billiard's position, velocity and velocity components at 
the collision points. 

The experiments of Feldt and Olafsen were motivated 
by the earlier work on the wedge billiard Q. Of course, 
that system is Hamiltonian so there is no dissipation or 
friction. Moreover, there is no ceiling or upper boundary, 
and the system is not driven. With these caveats in mind, 
it is instructive to compare the predictions of the original 
model with the actual experiments. These are displayed 
in Figures [5] and ini for a wedge with a half angle of 28.5°, 
exactly corresponding to the Feldt and Olafsen experi- 
ment. In Figure IS] we show a Poincare surface where the 
square of the normal velocity is plotted vs the tangential 
velocity after each boundary collision. As shown in the 
earlier work this choice generates an area preserv- 
ing map. The figure incorporates a number of distinct 
trajectories, all with a common energy. Surrounding a 
large, stable island associated with the period-one fixed 
point, we see a family of nested KAM tori all associated 
with the same fixed point. Surrounding this family are 
additional stability islands identified with different stable 
periodic points, as well as a space-filling chaotic orbit. 

In Figure [S] we display the same data plotted with the 
alternative coordinate pairs employed in the Feldt and 
Olafsen experiment. Clearly there is greater structural 
detail in the Hamiltonian version than in the experi- 
ment (Figure [5]) and simulations (Figure U). This is 
not surprising because, in the driven system, there is a 
distribution of energies. Moreover, the dominant role of 
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the fixed point is apparent in the former. In the plots 
of q2,n+ivs. 52, n for the simulations (Figure H]) the role 
of upper-boundary collisions is apparent. These are not 
visible in the published experimental work (Figure [S]) be- 
cause the scale height of the plots is too small. It is in- 
triguing that we see qualitatively similar behavior in the 
t„+ivstn plots for all three systems. Since the set of coor- 
dinate, q2/q2,mavVSU2/u2,max in the final picture closely 
corresponds to BirkhofF coordinates, the similarity with 
Figure [S] is not surprising. In future work it would be in- 
teresting to investigate the effect of the upper boundary 
on a Hamiltonian wedge. 
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FIG. 5: The spatial and temporal mappings of the collision 
heights (left column) and times of flight (right column) for 
the experimental data of Feldt and Olafsen From top 

to bottom: the parabola at 5.4 Hz, the wedge at 6.6 Hz, the 
hyperbola at 4.5 Hz the hyperbola at 5.8 Hz. 
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FIG. 6: A close-up of the spatial and temporal mappings of 
the collision heights (left column) and times of flight (right 
column) for the hyperbola (top row) and parabola (bottom 
row) at 4.5 Hz. 



FIG. 4: The spatial and temporal mappings of the collision 
heights (left column) and times of flight (right column). From 
top to bottom: the wedge at 6.6 Hz, the hyperbola at 4.5 Hz, 
the hyperbola at 5.8 Hz, the parabola at 4.5 Hz. 



VII. CONCLUSIONS 

There are open questions concerning how best to model 
impacts between systems of solid objects, such as granu- 
lar media. Examining the ergodic properties of a gravita- 
tional billiard provides an experimentally accessible sce- 
nario for testing and comparing a variety of impact mod- 
els. Here we have presented one model that captures the 
relevant dynamics required for describing the motion of 



a real world billiard for arbitrary boundaries. The model 
considers the more realistic situation of an inelastic, ro- 
tating, gravitational billiard in which there are retarding 
forces due to air resistance and friction. We have used the 
model to investigate driven parabolic, wedge and hyper- 
bolic billiards, and demonstrated that the parabola has 
stable, periodic motion, while the wedge and hyperbola 
(at high driving frequencies) appear chaotic. The hyper- 
bola, at low driving frequencies, has regular, periodic mo- 
tion, and behaved similarly to the parabola. The simple 
representation of the coefficient of restitution employed 
in the model resulted in good agreement with the recent 
experimental data of Feldt and Olafsen for all boundary 
shapes investigated iI3|, but not for secondary quantities 
derived from the data. The model also predicted addi- 
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FIG. 7: The normalized collision heights q2/q2,max versus 
the normalized tangential velocity Ui/u4,,max- Top row: the 
wedge at 6.6 Hz, the hyperbola at 5.8 Hz. Bottom row: the 
hyperbola at 4.5 Hz, the parabola at 4.5 Hz. 
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FIG. 8: Surface of section for the Hamiltonian gravitational 
billiard for a wedge half-angle equal to 28.5°. 



tional collisions not detected by the data. Gorski and 
Srokowski have also modeled the Feldt and Olafsen 
experiments. They employed a different model that in- 
cluded collisional energy loss, but ignored rotation. To 
achieve energy balance over long times and obtain qual- 
itative agreement with the experimental data, in their 
work it was necessary to invoke an unrealistic energy de- 
pendece of the coeficient of restitution. Perhaps this was 
due to their unconventional model which applied the re- 
duction in speed to the total velocity at collision, instead 
of the normal component (see Equation ([T2|) above) . The 
assignment of the value of the coefficient of restitution 
introduces the most uncertainty in modeling the billiard- 
boundary system, and resolution of this problem will re- 
quire additional experiments. It is interesting that the 
optimum numerical value is very different if rotation in- 
duced friction is not included. We will pursue this sur- 
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FIG. 9: Additional mappings of the Hamiltonian gravitational 
billiard for a wedge half-angle equal to 28.5°. The spatial and 
temporal mappings of the collision heights (top row), times of 
flight (middle row) and the normalized collision heights versus 
the normalized tangential velocity (bottom row) are shown. 



prising effect in a future work. 
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IX. APPENDIX 

Here we list the equations of motion for impacts of 
billiards with moving boundaries. The explanation that 
follows is a modified version of the impact theory set forth 
by Kane and Levinson for collisions between a sphere and 
stationary boundaries [18]. The theory considers the gen- 
eral three-dimensional case, but may be applied to two- 
dimensional systems by prescribing appropriate initial 
conditions for position and velocity. Since the numerical 
simulations presented in this paper are two-dimensional, 
we will point out simplifications when appropriate. 



A. Background 

Consider a sphere of mass m and radius b, whose mo- 
tion is confined by a moving boundary. The sphere has 
six degrees of freedom, three angles defining its orien- 
tation and three components defining its position. The 
boundary is infinitely massive and its shape is arbitrary. 
For the general billiard-boundary system, the inertial (or 
laboratory) frame is defined by three mutually perpendic- 
ular unit vectors (ni,n2,n3), where 112 is perpendicular 
to the plane formed by vectors rii and na, see Figure [TUl 
The reference frame at the collision point between the 
sphere and boundary is defined by the c— frame, and is 
related to the inertial coordinate frame by a translation 
and rotation of coordinates. It is more convenient to de- 
fine the collision response in a frame moving with the 
boundary (c- frame), oriented such that two of the unit 
vectors are locally parallel and orthogonal to the bound- 
ary surface at the collision point. Figure [TT] shows the 
c and n— frames for the two-dimensional case, where the 
frames are related by angle z. For curved boundaries an- 
gle z varies along the curve, and is uniquely determined 
for each collision. The sphere's angular and translational 
velocities at the collision point may be expressed in terms 
of the generalized speeds mi,...,M6, respectively, as 

u) = uiCi + U2C2 + U3C3 (15) 

and 

V = U4C1 -I- U5C2 + M6C3 (16) 

where v denotes the velocity of the sphere's center of 
mass. The boundary's velocity at the collision point is 
given by 

V = U4C1 + U^C2 + UqCs (17) 

where (for our simulations) U4, itg are found by taking 
components of v = ujAcos{ujt) ni. Note that for the 



two-dimensional case, the following generalized speeds 
are zero: ui, uq and Wg. 

From Equations ([T5l) and ([T6|. we define the velocity 
of the point P of the sphere that comes into contact with 
the boundary as 

= -v + u) X p 

= (U4C1 -I- U5C2 -I- M6C3) 

-|-(-UlCi -I- U2C2 + U3C3) X — 6C2 

— (U4 + 6U3)C1 + M5C2 + {uq — bui)c3 (18) 

Rewriting the result of Equation in terms of the 




FIG. 10; The inertial (or laboratory) coordinate frame. 
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FIG. 11: The coordinate frame at the collision point is de- 
noted by the c— frame. Reference frames c and n are related 
to each other by angle z. 
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generalized speeds, we have 

= (-6C3)U1 + (0)U2 + (6C1)U3 + (C1)U4 + (C2)U5 + (C3)U6 

(19) 

Then from Equation (|19p . the partial velocities of P, 
labeled as (r — 1,...,6), are determined by inspec- 
tion and are simply the coefficients for each generalized 
speed. For this problem ~ —bcs, = 0, Vg = 6ci, 

vj = Ci, = C2, = C3. 

Kane's impact model enables us to investigate the col- 
lision of a sphere as it impacts a boundary beginning at 
time ti and ending at time t2- Two dynamical equations 
essential to the model are the generalized impulse Ir and 
the generalized momentum p^. . The generalized impulse 
is generally applied to systems that are subjected to large 
action forces over a short time interval, and is defined as 



42 



R dt 



(r = l,...,6) 



(20) 



ti 



where v^(ti) is the partial velocity of the sphere at 
the point of contact with the boundary at time ii, and 
J*^ R dt is the contact force exerted on the sphere by the 
boundary at their contact point during the time interval 
[ti, t2\. Moreover, if we let S, ^ c; ■ J*^ R dt (i=l,2,3), 

then we can define J^^ R dt as 



t2 



Rdt = 5*101 + S2C2 + S3C3 



(21) 



ti 



The generalized momentum is defined as follows: 
dK 

Pr = ^ (r = l,...,6) (22) 



boundary. If there is no slipping at ^2, the following 
inequality must be satisfied: 



\t\ < /i |i/ 



(25) 



where T=SiCi + is the tangential impulse, i'=S'2C2 
is the normal impulse and fi is the coefficient of static 
friction. Consequently, 



C2 X (vs X C2) = 



(26) 



The equation states the tangential component of the ve- 
locity of separation is zero. If inequality ((25)) is violated, 
slipping occurs at t2, and t is expressed as 



C2 X (vs X C2) 

|C2 X (vs X €2)1 



(27) 



where the constant /i is the coefficient of kinetic friction. 

The model requires that the physical parameters b, m, 
J, e, /i and /i , and the generalized speeds at time ti are 
known, where J is the principal moment of inertia. Then 
the motion of the sphere at time t2 is fully defined by 
invoking Equations dig, ([25]), ^ and 



B. Connection Formulas 

The connection formulas define the relationship be- 
tween the billiard's pre-impact and post-impact veloci- 
ties for both the no-slip and slip cases for collisions on 
moving boundaries. For the no-slip case, the connection 
formulas are 



"2(^2) ~ U2(tl) 



(28) 



where K is the kinetic energy of the sphere and Ur are the 
generalized speeds. Integrating Equation (l20t results in 
the following approximation connecting the generalized 
impulse to the generalized momentum: 



' Pr{t2) -Pr{tl) 



(r = l,...,6) (23) 



U5{t2) = -eu^iti) + Ur^iti) [1 + e] 



(29) 



Ju^jti) + mb[u'^{ti) - M4(ti)] 
mo^ -I- J 



The approximation symbol in Equation (l23l) appears be- 
cause forces that remain constant during the time interval 
[ti, t2] are regarded as negligible. 

In order to capture the sphere's motion at time ^2, 
two assumptions supplement the use of Equation ([25)1 to- 
gether with a complete description of the sphere's motion 
at time ti. The first assumption is the normal compo- 
nents of the velocity of approach va and separation vs 
of the sphere, with respect to the boundary, have oppo- 
site directions, where the magnitudes are related by the 
following equation: 



C2 • vs = -eC2 • VA 



(24) 



In equation ([M)) . e is the coefficient of restitution. The 
second assumption determines if the sphere encounters 
no slipping or slipping at the point of contact with the 



^4(^2) = "4(^1) - buz{t2) 



(31) 



^ . Jui{ti) mh[uQ{ti) - Mg(ti)] 



ue{t2) = u^iti) + bui{t2) 



is! + s|y^'<^,\S2\ 



Si W m[u4{t2) - U4{ti)] 



S2 « rn[u5{t2) ~ uz{ti)] 



(33) 
(34) 
(35) 
(36) 
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(37) 



where e is the coefhcient of restitution, J is the princi- 
pal moment of inertia, m is the mass, b is the radius. Si 
are the impulses for i = 1,2,3 and fi is the coefficient 
of static friction. For no slipping successive use of Equa- 
tions (gll), dSOl), (ED), (ESI) and (ESI result in a set 
of values for ui,..., at time ^2- These values are valid 
if and only if inequality (|34|) is satisfied for values of , 
52 and ^3 given by Equations (05]), (06]) and dST]). For 
the two-dimensional case, wi, ug, Wg and 6*3 are equal to 
zero. 

Otherwise, if inequality is violated, then the 

sphere is slipping at time t2, and the quantities ^1(^2), 
^3(^2), ^4(^2), ""6(^2), 5*1 and S3 must be recalculated 
using the relationships listed below: 



7 = M6(il) - bui{tl) - Mg(tl) 



= i + ^ 

m J 



(38) 



(39) 



(40) 



Si « - 



A* " l'5'2| 



a|[l + (7/a)2]i/2 



5*3 « —51 



a 



"1(^2) ~ Mi(ii) - bSs/J 



U3{h)^U3iti) + bSi/J 



"4(^2) ~ '"4(ii) + Si/m 



ue{t2) w W6(ii) + Ss/m 



(41) 

(42) 

(43) 
(44) 
(45) 
(46) 



where a, 7 and fc are constants and /i is the coefficient of 
kinetic friction. Note that 1/2(^2), "5(^2) and S2 are given 
by Equations (EH]), (ESI) and respectively, regardless 
of whether or not the sphere experiences no slipping or 
slipping at time t2- For the two-dimensional case, iti, uq, 
Uq and S3 are equal to zero. 
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